/*
 *  process.h
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 15/11/07.
 *	Modified by Markus Lutz
 *	Modified by Sergio Conde on 17/06/2010
 *
 */

#ifndef _process_h_
#define _process_h_

#include "mecsystem.h"
#include "interface.h"
//#include "interfacearma.h"
//#include "interfaceeigen.h"
#include "ANALYSES/static.h"
#include <sys/stat.h>
#include "GEOMETRY/surface.h"
#include "history.h"
//#include <omp.h>


/*
#ifdef HAVE_SUPERLU
#include "sparsematrixslu.h"
#endif

#ifdef HAVE_UMFPACK
#include "sparsematrixumf.h"
#endif
*/

#include "MATH/sparsematrixslu.h"
#include "MATH/sparsematrixumf.h"
//#include "MATH/sparsematrixmumps.h"
#include <Eigen/Dense>

using namespace std;

extern int rankk;
extern int size;

class Interface;
//class InterfaceArma;
class InterfaceEigen;

class Process
{

public:
	int num;											/** Process number */
	int lrank;
	int dimension;										/** Space dimension */
	map<string, MecSystem*, less<string> > mMecSystems;	/** Map of systems belonging to this process. Key is name */
	MecSystem* rootSys;									/** Pointer to the root system of this process */
	map<string, int, less<string> > mSysProc;			/** Global table: system name - process number */
	vector<Process*> vPTProc;							/** Global Table: process number - pointer to process */
	bitset<20> relatedProc;								/** Bitset showing the relation of this process with others */
	vector<Interface> vInterfaces;						/** Interfaces of the process */
	//JCarma vector<InterfaceArma> vInterfaces;						/** Interfaces of the process */
	//JCeigen vector<InterfaceEigen> vInterfaces;						/** Interfaces of the process */
	multimap<string, int, less<string> > localSharedPoints;	/** Local Point name - host process */
//	Analysis* lastAnalysis;								/** Last performed analysis */
//	Analysis* currentAnalysis;							/** Current analysis being performed */
	int numDofs;										/** Number of dofs of the process */
	vector<double> q;									/** position dofs */
	vector<double> dq;								/** velocity dofs */
	vector<double> ddq;							/** acceleration dofs */ // edit
	vector<double> res;							/** residual vector */
	vector<double> qConv;							/** converged q */
	vector<double> dqConv;							/** converged velocity */
	vector<double> Deltaq;							/** increment in position dofs, from the last converged value */
	bool isVelNeededLast;								/** 0 if velocity was not needed at the last analysis */
	map<string, Material*, less<string> > mMaterials;	/** Map of materials */
	map<string, Surface*, less<string> > mSurfaces;	/** Map of surfaces */
	char outFileName[30];
	char title[100];
	ofstream ou;
	vector<int> localGlobI;								/** localGlob for interior dofs */
	vector<int> globalLocalI;							/** Input: global dof / Output: interior dof (-1 if at the boundary)*/
	vector<int> localGlobB;								/** localGlob for boundary dofs */
	vector<int> globalLocalB;							/** Input: global dof / Output: boundary dof (-1 if at the interior)*/
	int numDofsI;										/** Number of interior dofs */
	int numDofsB;										/** Number of boundary nodes */
	sparseMatrix* Kint;							/** Tangent of interior dofs */
	vector<int> vDofsOut;								/** Global dofs numbers of points with requested output */
	vector<int> vStressOut;								/** Global dofs numbers of points with requested output */
	char outDofsFile[30];
	char outVelDofsFile[30];
	char outGraphFile[300];
	char outTempFile[30];
	char outIterFile[30];
	char outStressFile[30];
	ofstream fDofs;										/** Output stream for dofs of requested points */
	ofstream fVelDofs;										/** Output stream for dofs of requested points */
	ofstream fGraph;									/** Output stream for graphical output */
	ofstream fTemp;										/** Output stream for temperature of requested bodies */
	ofstream fIter;
	ofstream fStress;
	int vecType;										/** LMX vector type */
	int matType;										/** LMX matrix type */
	string procSS;			/** string containing process-numStep */
	double grav;	// Gravity value
	int gravDir;	// Gravity direction
	double eneResFi;
	double eneResFq;
	double eneResFt;
	double eneResInter;
	double eneResMax;
	double eneRes;	
	vector<double> resFq;
	vector<double> resFt;
	vector<double> resInter;
	vector<int> sparseIntTot;	// relates positions in a vector for sparse format
	sparseMatrix* M;
	bool dynFlag;
	double Ekin, Epot, Econ, Eamb, Eext, Dint;
	double Entropy, EntropyAmb;
	vector<double> AM;
	vector<double> vCon;
	int numDofsInertia;	// number of dofs with inertia
	vector<int> numNoInertia;	// global dof numbers that do NOT contribute to inertia
	double thetab;	// Ambient temperature
	vector<string> vTempOut;	/** Names of bodies/joints with requested output */
	bool projFlag;
	int autoProjFlag;
	//	Juan Luis
	sparseMatrix* MfacExp;
	vector<double> vLumpedMass;
	sparseMatrix* MfacProj;
	double projPenaltyMax;
	vector<int> ia;
	vector<int> ja;
	vector<int> iaInt;		// Row index (internal tangent matrix)
	vector<int> jaInt;		// Column index (internal tangent matrix)		
	vector<int> iaInertia;			// Row index (global mass matrix)
	vector<int> jaInertia;			// Column index	(global mass matrix)
	sparseMatrix* pMatrix;	// P = (A + alpha*Phiq*Phiq) such that P*qdot = A*qdot_ori
	vector<double> projVel;
	vector<vector<double> > Phiq;
	// Projection stuff
	vector<vector<double> > Gamma;
	vector<vector<double> > GammaFac;
	vector<vector<double> > auxGamma;
	//JC
	vector<vector<double> > Kbb;	// Boundary dofs tangent
	vector<vector<double> > Kib;
	vector<vector<double> > B;
	vector<vector<double> > Kbi;
	vector<double> Ri;
	vector<double> Rb;
	bool symmFlag;
	multimap<string, string, less<string> >	mPointsNoDofs;
	vector<double> vPointsEndTimes;
	vector<int> vNoDofs;	// global dof numbers to remove because are fixed
	vector<double> vNoDofsEndTimes;	
	sparseMatrix* MSmooth;
	vector<vector<double> > FSmooth;
	vector<vector<double> > FSmoothDerived;
	MatrixXd SigmaMat;
	vector<double> PSigma;
	double VMStress;
	multimap<string, vector<string>, less<string> > mPointsImposedDofs;
	string dispFile;
	vector<int> vImposedDofs;		// imposed global dof numbers
	vector<string> vImposedDofsHis;	// names of histories, corresponding to the imposed dofs
	//map<int, string, less<int> > mImposedDofs;
	map<string,History*,less<string> > mHistories;
	vector<double> qIni;	// Initial dof vector
	char timeInfoFileName[30];
	ofstream ouTime;
	double tini;
	int numInterfaces;
	SuperMatrix	B_temp;
	SuperMatrix	X_temp;
	
public:
	Process();
	Process(int, int, int, int, int&, const char *);
	~Process();
	
	int readSystem(string, ifstream &, int);
	void setSysProc(map<string, int, less<string> >, vector<Process*>);		/** set global tables that relates systems-processes */
	Process* getPTProc(string);							/** Gets a pointer to the process that contains an specific system */
	Process* getPTProc(int);							/** Gets a pointer to the process number num */
	MecSystem* getPTSys(Process*, string);				/** Gets a pointer to a system, from a pointer to its process */
	MecSystem* getPTSystem(string);						/** Gets a pointer to a specific system, that may belong to a different process */
	void buildSysTree();								/** System tree of the process (not global) */
	void createPoints();								/** Create points within the systems */
	void setPPointersCompBodies();						/** Set pointers to points at body components */
	void setPPointersLoadBodies();						/** Set pointers to points at load components in bodies */
//	void setPPointersCompJoints();						/** Set pointers to points at joints */
//	void setPPointersLoadJoints();						/** Set pointers to points at load components in joints */
	void setPPointersConsJoints();						/** Set pointers to points at constraint in joints */
	void setPPointSurface();							/** Set pointers to surface points in bodies */
	void setRelProc();									/** Set crossed relations among processes */
	vector<Point*> getPTPoints(vector<string>&);		/** Gets local pointers from a vector of point names */
	Point* getPTPoint(const string&);					/** Gets local pointer from a single point name */
	void createInterfaces(const int&);					/** Create interfaces at this process with other processes */
	void relateInterfaces();							/** Relate interfaces*/
	void calcSizeInterfaces(const int&, bool&);		/** Calc proper dof size to interfaces */
	void setSizeInterfaces(bool&);						/** Set proper dofs to interfaces */
	void addRemoteSharedPoint(string, int);					/** Adds to the remote process the name of the point and its process number */
	void addLocalSharedPoint(string, int);					/** Adds to localSharedPoints the name of the point and the host process */
	void setDofsBodies();									/** set dofs at points of bodies */
	void setDofsRoots();									/** set dofs at points at root systems */
	void setDofsJoints();									/** set dofs at points at joints */
	void setGlobalDofNumbers();								/** setup the global dof numbering of this process */
	void setLocalGlob();									/** Set local-global numbering on components */
	void sendLoadToBody(string, string, string, ifstream &, const vector<double> &);	// edit /** Send load to body */
	void initializeGlobalDofsBodies();							/** Initialize dofs at bodies*/
	void synchronizeSharedDofs();							/** Initialize dofs that have not been initialized by any component (joints) */
	void initializeGlobalDofsJoints();							/** Initialize dofs at joints*/
	void setVectorQ(vector<double>& q_in);				/** Set vector q=q_in */
//	void checkVelNeeded();									/** Check if velocity vector is needed */
	void setMaterials(const map<string, Material*, less<string> >&);	/** Set materials at the process */
	void setMatAtComp();								/** Set pointers to material at components */	
	void formIsolatedResidue(const double& t, const int&, const int&, const int&, const double&);			/** Form residue from the isolated process */
	void formIsolatedTangent(double& t, const int&);			/** Form tangent from the isolated process */
	void formIsolatedTangent(double& t, double&, double&, const int&, const int&, const double&);	/** Form tangent from the isolated process */ // edit
	void resetInterInt();									/** Reset internal interface contributions */
	void factorizeInternalMatrix(int);							/** Factorize tangent of interior dofs */
	void setInterExt(bool);										/** Set external interface contributions */
	void setInterInt();										/** Set internal interface contributions */
	void formResidue(const int&);							/** Form residue, with contributions from interfaces */
	void formTangent();							/** Form tangent, with contributions from interfaces */
	
	void myEvaluationImplicit	(const double& time, const int&, const bool&);		/** Evaluation for implicit integrators */
	void myEvaluationExplicit	(const double& time, const bool&, const bool&);		/** Evaluation for explicit integrators */
	/*
	void myEvaluationExplicit	(	const vector<double>& q,
						const vector<double>& qdot,
						vector<double>& qddot,
						const double& time
						);
	*/
	void solveUpdate(vector<double>&,int, double&);	//Solve tang*x = res and update

	void reportContent();
	void reportPointsDofs();
	void reportPointsDofsFile();
	
	bool checkIfSystem(string sysName)	{ if(mMecSystems.count(sysName)) return(true); else return(false); }/** Check if a system is within the processs */
	int getNumberInterfaces();
	void createDofOutFile();				/** Create an output file to store requested dofs */
	void createStressOutFile();				/** Create an output file to store requested dofs */
	void fillVDofsOut(const string&);		/** Fill the vector vDofsOut, given a local point name */
	void fillVStressOut(const string&);		/** Fill the vector vDofsOut, given a local point name */
	void writeFileDofs(const double&);
	void writeFileStress(const double&);

        void writeFileGraph(string &, string &, const int&, bool&, bool&, bool&);
	
	void writeToFile(const string &);
	void writeToFile0(const string &);
	void writeToFile(const string &, const char *);
	void writeToFile(const string &, const string &, const string &);
	void writeToFile(const string &, const int&);
	void writeToFile(const string &, const int&, const string &, const double &);
	void writeToFile(const int &);
	void writeToFile(const int &, const string &);
	void writeToFile(const vector<double>);
//	void writeToFile(const lmx::Matrix<double>);
	void writeToFile(const string&, const int&, const int&, const string&);
	void writeToFile(const vector< vector<double> >&);
	void writeToFile(const vector<int>&);
	void writeToFileTime(const string& , const double&);
	void setDofValues(vector<int>&, vector<double>&);	//Set dof values at given positions
	void initializeRes(const int&, const bool);				/** Initialize residual vectors.  st=0 Static, st=1 dynamic*/
	void initializeTang(const int&, const bool, const int, const bool, const bool);				/** Initialize tangent matrices. st=0 Static, st=1 dynamic*/
	void setSparsePattern();	/** Computes the sparse pattern of global and internal tangent matrices */
	void setSparsePatternInertia();	/** Computes the sparse pattern of global and internal tangent matrices */
	void setGravity(double, int);

        void createGraphOutFile(const string &, const string &, const int&);			/** Create an output file to store graphical output */
	void closeGraphOutFile(const string &, const int&);				/** close output graphic file with specific commands (needed by pvd format) */
	void closeFileDofs(void);
	void closeFileStress(void);
	void setQConv(void);		// Store the converged step, used in dynamics
	void calcEner(void);
	void sendIniPosToPoint(string&, string&, vector<double>&);
	void sendIniPosToPointSet(string&, string&, vector<double>&);
	void sendIniVelToPoint(string&, string&, vector<double>&);
	void sendIniVelToPointSet(string&, string&, vector<double>&);
	void sendIniVelToBody(string&, string&, string&, vector<double>&, vector<double>&);
	void setInitialConditions(double&);
	void update(const int&);
	void setInertiaDofNumbers();	// Dofs associated with inertia (they contribute to mass matrix)
	void createConvergedVectors();
	void removeFromVector(vector<double>&, vector<int>&, vector<double>&);
	void restoreFromVector(vector<double>&, vector<int>&, vector<double>&);
	void setTemperature(void);
	void createTempOutFile();				/** Create an output file to store requested temperatures */
	void fillVTempOut(const string&);		/** Fill the vector vDofsOut, given a local point name */
	void writeFileTemp(const double&);
	void closeFileTemp(void);
	void sendIniPosToBody(string&, string&, vector<double>&);
	void createIterOutFile();	// Creation of file to store iterations
	void writeFileIter(const double&, const int&);
	void closeFileIter(void);
	void calcEntropy(void);
	void calcAM(void);
	void calcCon(double&);
	void setQatComponents();
	void updateLagMul();
	void checkConvLagMul(bool&);
	void projectVel(double&, int);
	void setPatternInertia(sparseMatrix*);
	double calcProjPenalty(double&, string&);
	double calcDeltaEpRes(double, vector<vector<double> >&, vector<double>&, double);
	double calcDotPhiRes(double, vector<vector<double> >&, vector<double>&, double);
	void setSurfaces(const map<string, Surface*, less<string> >&);	// Set surfaces at the process
	void setSurfAtBodies();								// Set pointers to surfaces at bodies
	void createContacts();
	void addNoDofsPoint(string&, string&, double&);
	void setNoDofs();
	void createLumpedMass();
	void synchronizeLumpedMass();
	void repairChangedPointNames();
	void fillMSmooth(double &);
	void calcSmoothStress(double &);
	void addImposedDispPoint(string&,string&,char*);
	void addImposedDispPointSet(string&,string&,char*);
	void setImposedDofs();
	void setHistories(map<string,History*,less<string> >&);
	void setQIni();

	int checkComponents();
	int setComponents();
	
	// Friend functions
//	friend void lu(lmx::Matrix<double>&);
	friend void printSTLVector(vector<int>& v);
	friend void printSTLMatrix(const vector< vector<int> >& m);
	friend void printSTLMatrix(const vector< vector<double> >& m);
	friend void printSTLVector(vector<double>& v);
//	friend void forSub(lmx::Matrix<double>& mat, lmx::Vector<double>& vec);
//	friend void backSub(lmx::Matrix<double>& mat, lmx::Vector<double>& vec);
//	friend void forSub( lmx::Matrix<double>& mat, lmx::Matrix<double>& mat_r);
//	friend void backSub(lmx::Matrix<double>& mat, lmx::Matrix<double>& mat_r);
	friend void setZero(vector<double>& v);
	friend void setZero(vector< vector<double> >& m);
	friend void printSparseMatrix(sparseMatrix& m);
	friend void lu(vector<vector<double> >&);
	friend void forSub(vector<vector<double> >& mat, vector<double>&);
	friend void backSub(vector<vector<double> >& mat, vector<double>&);
	friend SuperMatrix convert_stl_SLU(vector<vector<double> >&);
	friend void convert_SLU_stl(SuperMatrix&, vector<vector<double> >&);

};

#endif
